#include "func_secondderr.h"
#include <malloc.h>


struct iteration {
	boolean active; // If the record holds any data
	uint32_t N;
	double max;// = max |d|
};


#define RESCALING_2_FOR_OK 5 // If result is the same within `eps` after multilying knot number N by 2^{this}, that's SUCCESS

int func_secondderr_max(complex (*func)(double , void *), void * func_context, double a, double b, double eps, double * result, uint32_t * nknot){
	uint32_t N = 300; /* Knots at [a,b] including borders */
	int nitr = RESCALING_2_FOR_OK + 1;
	struct iteration itrs[nitr];
	{
		int j;
		for(j = 0;j < nitr;j += 1){
			itrs[j].active = FALSE;
		}
	}
	for( ; ; N = N * 2){
		//printf("%i\n", N);
		if( N >= UINT32_MAX / 2 ){ //Futuristic limit
			return 2;
		}
		struct iteration new_itr;
		new_itr.active = TRUE;
		new_itr.N = N;
		
		uint32_t o;
		double h = (b - a)/(N - 1);
		// Calculating maximum of second derrivative absolute value
		double max = 0;
		complex prev[2];//@ new iteration, previous values: 0 - previous, 1 - preprevious
						//uninit until o = 2, of couse
		prev[0]=0;
		prev[1]=0;//Just in case
		for(o=0;o<N;o+=1){//Where o is knot id from 0 to N-1
			double x = a + (o*(b-a))/(N-1); 
			complex here = (*func)(x,func_context);
			if(o >= 2) { // calcule second derrivative mod and make comparations
				complex derr2 = (here - 2*prev[0] + prev[1]) / (h*h);
				double here_mod = cabs(derr2);
				if(here_mod > max){
					max = here_mod;
				}
			}
			//push
			prev[1]=prev[0];
			prev[0]=here;
		}
		// Saving 
		new_itr.max = max;
		// Shifting iterations
		int ksi;
		for(ksi = nitr; ksi > 1; ksi -= 1){
			itrs [ksi - 1] = itrs [ ksi - 2 ];
		}
		itrs [0] = new_itr;
		
		// Checking success condition
		// We have success, if the difference of `max` between this 
		// iteration and all the previous is within specified `eps`.
		// This means, the difference has been no more than 2*eps
		// between every value while shrinking knot step from initial
		// to 2^(RESCALING_2_FOR_OK).
		//
		// It seems that (a-2b+c)/h^2 method is O(h^2)
		// N ~ 2^itr, h~2^(-itr), diff ~ 1/2^2*itr, 2^(-2*RESCALING_2_FOR_OK) << 1
		// ~ the method seems to realy estimate if we are close enough
		// to true answer
		boolean precision_achieved = TRUE;
		for (ksi = 1; ksi < nitr - 1; ksi+=1){
			if ( ! itrs[ksi].active ){
				// Full test not yet ready 
				precision_achieved = FALSE;
				break;
			}
			double diff = fabs(itrs[ksi].max - itrs[0].max);
			double rdiff = diff/itrs[0].max;
			if ( rdiff < eps ){
				continue;
			}else{
				precision_achieved = FALSE;
				break;
			}
			// Zero f`` function will nerer trigger break condition, the program would consume memory and fail
		}
		
		if ( precision_achieved ){
			* result = itrs[0].max;
			break;
		}
	}
	if ( nknot != NULL) {
		* nknot = N;
	}
	return 0; /* Success */
}
